Many-body correlations in Semiclassical Molecular 
Dynamics and Skyrme forces for symmetric Nuclear 
2 : Matter 



o 

+-» 

O 

O 

(N 



(N 



M Papa 

INFN - Sezione di Catania, Via S. Sofia 64, 95123 Catania, Italy 
E-mail: massimo.papa@ct.infn.it 



Abstract. Constraint Molecular dynamics CoMD calculations have been performed for 
symmetric nuclear matter (NM) by using a simple effective interactions of the Skyrme type. The 
set of parameter values reproducing common accepted saturation properties of nuclear matter 
have been obtained for different degree of stiffness characterizing the iso-vectorial potential 
• density dependence. A comparison with results obtained in the limit of the Semi- Classical 

| Mean Field approximation performed using the same kind of interaction put in evidence the 

£^ ■ role played by the many-body correlations present in the model explaining also the noticeable 

differences obtained in the parameter values in the two cases. 
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1. Introduction 

. The description of many-body systems is one of the most difficult problems in nuclear physics due 

to the complexity of this kind of systems which are quantum objects described by a large number 
of degrees of freedom. A large variety of theoretical models have been developed using mean-field 
based and beyond-mean field approaches like Density Functional Theory pQ and Energy Density 
Function Theory [2J [3] . In these approaches which use the independent particle approximation as 
a starting point, phenomenological effective interaction like Skyrme and Gogny forces are widely 
used taking advantage of their simple form[5|. In nuclear structure modeling we can quote the 
Skyrme-Hartree-Fock (SHF) methods, the relativistic mean-field (RMF) approach (as well as 
their Bogoliubov extensions), and the Hartree-Fock-Bogoliubov method with the finite-range 
Gogny force [5j [6] . 

Complexity still becomes higher when nuclear dynamics, triggered by nuclear collisions, is 
studied to understand the property of nuclear forces far from the stability. At the Fermi energy 
and beyond semiclassical methods become necessary to describe the produced processes in which 
practically all the degrees of freedom are involved. Many-body correlations are responsible for 
the reorganization of the hot Nuclear Matter (NM) in to clusters. Two main classes of approaches 
have been developed to handle this complex scenario. The first one is based on the Boltzman 
transport equation including, at the higher energy, also three body collisions term [HIE]- In this 
case the theoretical approach is able to describe the time evolution of the one-body distribution 
function in phase-space. Apart from the collision term, the Boltzman equation corresponds 
to the semi-classical limit of time-dependent Hartree-Fock equations in phase-space. Different 
models have been implemented on this ground. They essentially differ from each other for the 



strategy adopted to simulate the collision terms and including the method used to represent the 
phase-space distribution starting from the nucleonic degree of freedom (test particle methods). 
As in Harthrre-Fock mean field methods the main ingredient concerning the interaction is the 
energy density which for a Skyrme type phenomenological interaction has a rather manageable 
form. In particular the local part of the interaction produces a simple functional for the potential 
energy expressed through an algebraic function of the density (see also Sec. II). In the early 90 
a further development included stochastic forces (typical of Langevin processes) |10] in the so 
called Stochastic Mean Field[TT] model to produce fluctuations capable of describing associated 
phenomena like mean-field instabilities leading to the cluster production. The second line of 
development of theoretical semiclassical approaches to deal the nuclear many-body problem is 
represented by the so called Molecular Dynamics models. The point of view of these methods 
is some sense antithetic compared to the ones based on the mean-field concept. However also 
in these cases, due to their simplicity the same kind of phenomenological effective interactions 
Skyrme and Gogny are widely used. The starting point of all these models is a basic assumption 
on the wave function describing the single nucleon which is represented through a wave packet 
well localized in phase-space with uncertainty satisfying the uncertainty principle. The many- 
body wave fucntion can be expressed through a direct product, leading to the so called Quantum 
Molecular Dynamics (QMD)-like models \12\ fl"3j [Tl"l [T5] or an anti-symmetrized product like 
in FMD and AMD |16t I17j . Variational principles have been used to obtain the equation of 
motion for the the many-body functions. Within their own approximation schemes these models 
treat the many-body problem without the usage of the mean-field approximation. Many-body 
correlations are spontaneously produced and in general are able to describe the main features 
concerning the multi-breakup processes [18] observed in heavy ion collisions at intermediate 
energy and beyond. Due to these deep conceptual differences between the mean-field based 
models and the molecular dynamics ones it becomes interesting, and necessary in our opinion, 
to investigate the differences between the energy-density functionals which are produced with 
these two classes of approaches when one uses, in both the cases, the same kind of effective 
elementary interaction between nucleons. The study presented in this work is performed by 
using the Constraint Molecular Dynamical Model CoMD [13, 14] . It is a part of a forthcoming 
more extended study involving also asymmetric NM. As it will be shown the presence of iso- 
vectorial forces strongly affects the results of this comparison (see also [E]). The quantitative 
results obviously will be related to some specificity of the model and to the used effective 
interaction, but as due to the common features which links most of the molecular dynamics 
models (semi classical wave packet dynamics), the obtained results could assume a more wide 
meaning at a qualitative level. The work is organized as follows: in Sec. II we illustrate the 
choice of the effective interaction and the related total energy functional are introduced. In 
Sec. Ill we describe the NM simulations. The results are presented in Sec. IV. Sec.V contains 
the summary and the concluding remarks. 

2. The effective interaction and the total energy 

Before illustrating with some detail the model and the NM simulations, in this section we 
briefly introduce the effective interaction which we use in our calculations. We also present the 
way in which the related total energy is obtained in the case of the semiclassical mean-field 
approximation (Se-MFA) and in the case of the semiclassical-wave packet dynamics. 

The elementary interaction between two nucleons with spacial coordinates r,r' and third 
isospin component r,r' is of the Skyrme type and has the following form: 



V(r,r') = V^5(r-r') 



Zk*( r _ r ') + 5{r 



r>) + ^F>(26 T y-l)6(r-r>) (1) 
Po 



The first term of eq.(l) represents the iso-scalar contribution to the two-body interaction, the 
second one is the usual 3-body effective interaction. For this term the p density dependence 
is modeled through the parameter a. This term appears to be a generalization of the Brink- 
Voutherin [3j three-body term corresponding to a = 2. This generalization is sustained both on 
the basis of phenomenological parametrization used in mean-field approaches [8J and on more 
complex microscopic calculations based on G matrix calculations with realistic interactions. 
The third term represents instead the Iso- Vectorial contribution. The form factors Ff. have the 
following expression: 

F k = {p/ PQ )F' k (2) 
, 2(p/ Po ) 

Fl ~ TTp7p- (3) 

F' 2 = 1 (4) 
K = (P/P0)- 1/2 (5) 

n = (p/poV- 1 (6) 



po is the saturation density value. Also in this case the Fk form factors are introduced to 
reproduce the results of more complex microscopic calculations based on realistic interaction 
\20\ [2T| [22j [23"1 [23] concerning the symmetry energy and reflecting effects beyond the two- 
body interaction . These form factors have been widely used in Hartree-Fock calculations, in 
semiclassical BUU calculations and in others molecular dynamcs models. In particular, F4 will 
be used, in the limit p/po ~ 1, to perform calculation with stiffness parameter 7 different from 
the ones related to the others form factors (in the same limit F\ , Fi , F3 correspond to 7 values 
1.5, 1 and 0.5 respectively). Finally we observe that for simplicity reason we do not add to this 
interaction non local effects. Regarding this assumption we note that some of the suggested 
Skyrme parametrizations have a vanishing or very small terms related to the non locality at the 
saturation density in the limit of Nuclear Matter (see for example Ref. [25] ) . Moreover as it will 
be shown in the following, in CoMD calculations the effective potential which determines the 
effective force able to govern the motion of the wave packet centroid is already non-local (sum 
of Gaussian contributions depending on the intra nucleons distances) being the results of the 
convolution of the well localized wave packets with 5 function in space. 

Even if well-known, we shortly recall in the next section the main ingredients of the Se-MFA. We 
think this section useful to present in the following the way in which this assumption is instead 
broken in the semiclassical wave-packets dynamics. 



2.1. The Semiclassical Mean Field Approximation 

At a fixed time starting from eq.(l) (we are considering here a stationary problem) , in a 
rather general way, we can obtain the expression for the total energy W related to the potential 
interaction by folding the elementary interaction with the two-body density distribution in 
phase-space D(r, r' , p, p'). By taking into account the usual truncation condition on D typical 
of the mean-field approximation, we can write: 

D = Z>!(r,p) ■ Z^i/, pO (7) 

where D\ is the one-body distribution. Therefore by using eq.(l) for W we get: 

W = - J V(r,r')D(r,r',p,p')drdr'dpdp' = -( V (2) D 1 (r,p)D 1 (r,p')drdpdp' (8) 



Taking in to account that by definition / Lq(r,p)dp = p(r) we obtain: 



W = i /" V^p 2 dr (9) 

appear to be the energy density associated to the potential energy from which we can 
obtain the related binding energy per nucleon due to the potential interaction E pot 

E P ot = U twb + U trb + U asy = \v^p =\—+ . T f° - + ]-T 4 F k (p)p 2 (10) 

I 1 p [a + i)p z 

f3 = pn p Pv represents the charge/mass asymmetry parameter evaluated from the neutron and 
proton density p n , p p respectively. The total binding energy E is obtained by adding the kinetic 
contribution coming from the Fermi motion 

E = Epot + E^in (11) 

In E^ in j3 terms of order greater than two are neglected. In particular we see how the iso- 
vectorial vectorial forces with strength proportional to T4 contribute to the symmetry energy 
E sym depending on the asymmetry parameter (3. For quadratic form in /3 we get: 

Esym{p) = dsymP (13) 



1 d 2 E 1 h 2 3n 2 p 2/3 1 

^sym- 2 ^ d(3 2)- 6m ^ 2 ) + 2 



( = 'a— (^) 2/3 + n^Fip) (14) 



The common accepted bulk value of e sym at the saturation density is about 30 MeV even if 
relativistic Hartree-Fock models can predict higher values up to about 40 MeV |26tl27|. However, 
still under study is the density dependence of this quantity which is able to affect neutron-skin 
thickness in nuclei and the Giant Monopole Resonances |26t |2"71 12B] . Finally we note that the 
structure of the obtained energy-density functional in Se-MFA makes independent the choice of 
the parameter describing the main properties of the symmetry energy from the other ones which 
instead are fixed from the saturation properties of the symmetric nuclear matter. We will show 
in the next section that this is no longer true in the CoMD approach. 

In molecular dynamics approaches a strong assumption is done on the wave functions 
describing the nucleonic degree of freedom. It is commonly assumed that the wave function 
is a gaussian wave packet with fixed width a r in coordinate space. The centroid in phase-space 
is indicated with r$, Pi 

T ex Pi 7T3 ( 15 ) 



The Wigner transform of (pi is 



r^I-^-^l (16) 



the widths in momentum and space satisfy the minimum uncertainty principle condition 
a r a p = hh. Another assumption concerns the N-body Wigner distribution that in CoMD model 
|13J (like in QMD approach [12]) is a direct product of the single particle distributions. In this 



case therefore for a system formed by A nucleons the one-body and 2-body distributions above 
introduced have a multi-component structure that is: 



A 

L>i(r,p) = £/i(r,p) (17) 
1 

A 

D(r,p,r',p') = J2 Mr,p)f,(r',p') (18) 

From the above relations we can see that in general D ^ Di(r, p)_Di(r', p'). In particular, for 
the case in which these distributions are expressed as a sum of different localized components 
inside a volume V g , we can indicate with a g the number of components which give non negligible 

contributions in V g . The relative difference associated to eq.(18) 1 — £>1 ^ r ' p ^ 1 ^ r ' p - can be 
estimated to be of the order of l/a g which corresponds to the ratio between diagonal (i = j) 
and off-diagonal elements (i ^ j) within the ensemble of a g (a g — 1) terms. In semiclassical mean- 
field models a g can be enough high to make the difference negligible, in fact the single particle 
distribution usually spreads over the whole system (test particles methods) and therefore the 
truncation condition eq.(18) can be retained valid El |9]. On the contrary this is not surely 
the case for the molecular dynamics approaches for which the typical spreading volume is of the 
order of 2-10 /m 3 . Localization and therefore coherence of the wave-packed used to describe 
the single-particle wave-functions allows to keep memory of the two-body nature of the inter- 
particles interaction, and at same time, allows for the spontaneous appearance of the clustering 
phenomena in simulations concerning low-density and excited portion of nuclear matter. With 
these assumptions on the 2-body phase-space distribution, taking in to account the properties 
of the 5 function, we can obtain the explicity expression for the different terms concerning the 
total energy; for the two-body isoscalar contribution Wt w b we get: 

w** = o r/ 2W2) E «q»[- ^fffi ( 19 ) 

2^(47^2)3/2) 4cj2 

W twb = t^£S* (20) 

In the above expression S l v is the normalized sum of the Gaussian terms and it represents 
just a measure of the overlap between the nucleonic wave-packets. Its two body character is 
quite explicit. In the calculations concerning the NM simulation that we will illustrate in the 
next sections , the large number of particles A involved in the system allows us to write the 
above quantity in a simpler way by introdcing the average overlap per nucleon Si = S v whose 
dependence on the particle index in the ideal case can be omitted: 

W twb = ^S~ v (22) 
„ _ W tw b _ T — 

&twb — -. — — ~ O v [Z6 

A 2p 

By comparing the expression obtained for Et wb in the two different appraches (eq.(10) and 
eq.(23)) we note a formal analogy where the variable p is substituted by the overlapp integral S v 



per nucleoli. We however observe that this analogy is only formal, this aspect will be discussed in 
some detail in the next subsection. Concerning the three-body term, according to the evaluations 
reported in reference [12J and taking in to account the previous obervations we get: 

Etrb = (a+lK^ (24) 
For the term related to the iso-vectorial interaction and for the most simple case F' = 1 in the 
limit A,N, Z >> 1 (iV and Z represent the number of neutron and protons) we obtain: 

W isv = I^_( N 2pnn + Z 2^pp _ 2NZ ~ p np^ ^ 
2/50 

where p cc with cd equal to nn, pp and np represents the overlap integral per couples of neutrons, 
protons and neutron-proton. A more convenient form for the above expression is obtained by 
introducing the two followng quantities: 

_ + Z^v 

9 ~ N* + Z* [ZU) 



a = ' - ' (30) 
P 



for symmetric NM N = Z and we get: 

Wg = -^A 2 pa (31) 

E?sv = -^-F^)p A a = E Mas (32) 
4po 

where pa = Ap. The expression in eq.(32) also contains a generalization to the cases in which we 
use the generic form factors Fj.. Here F^ keep the same functional form in S v using the formal 

analogy ^ — > == above discussed. From eq.(32) we obtain that for a ^ (as CoMD calculations 

predicts, see next sections) the iso-vectorial force produces a term which we name Eyi as . The 
related effect is not-negligible if compared to the balance of the different terms appearing in the 
expression of the total energy. The correlation coefficient a by definition (see eq.(30)) represents 
the difference in percentage of the overlap between the neutron-proton couples from the one 
related to the couples formed by homonym nucleons. It also depends on the strength of the 
iso-vectorial forces (T4 parameter). 

Finally, the kinetic contribution E^in is obtained in a self-consistent way for ground state 
configurations by applying the cooling-warming procedure coupled with the constraint on the 
occupation numbers (see Ref.}14j) given by the Pauli principle. The total energy per nucleon 
is therefore obtained by adding all the contributions discussed in this section and it will be 
indicated as E c = E, pot (S v , a, pa) + Eki n {p)- In particular we note that, at this level, it depends 
on the new defined primary quantity S v , a, p~A- In the following section we will try to relate 
these quantities to more fundamental ones that are the density p and the spatial correlation 
function v between nucleon pairs. 



2.2. Overlap integrals and spacial correlations 

In a rather general way and as suggested from CoMD calculations (see Fig.l), for an uniform 
many-body system at a fixed density we can introduce the probability p to have a particle in 
the volume dV\ localized in ri and a second one in the volume dVi localized in x<i- Due to the 
uniformity condition p depends only on r = |i"i — r2|. p can be expressed in the following form: 



with i/(0) = 1 and ko > 0. Moreover lim r _ 5>00 z/(r) = 0, that is: no spacial correlations can be 
expected for very distant particles. Non zero values of v can be instead expected at relatively 
small distances due both to the interaction (for an attractive interaction the positive sign must 
be considered in eq.(33)) and to the symmetry of the many-body wave function for quantal 
systems of identical particles (in the case of identical Fermions we must consider the sign minus 
and ko = 1 ) . For a classical system of non interacting particles we have a vanishing correlation 
effect p = 1 (ko = 0). For our aims, we need to evaluate a normalized probability P = cp in 
such a way: 



Vc represents the volume in which the well localized correlation function v is different from zero. 
This volume is always finite and of the order of the a r 3 in our model calculations. So that in 
the limit V — > oo we get c = p and : 



In Fig.l, just as an example, we display results of calculations for a given set of parameters for 
the Skyrme interaction (see the next section) . The calculations show the probability to find two 
nucleons at a distance r in case of Pauli blocked couples Pi (red points, neutron and proton 
couples with same spins), for the case of un-blocked proton and neutron couples Po (black point) 
and finally for neutron-proton couples P np (blue points). 

The calculations are referred to a symmetric portion of nuclear matter consisting of a sphere 
containing 2000 nucleons at a density of about 0.17/m -3 . The calculations include also the 
iso- vectorial potential energy (T4 = 59 MeV in this example ). The probability distribution has 
been estimated by taking in to account only couples of nucleons in which at least one of the two 
nucleons is localized inside a smaller sphere, centered at the origin, with a radius of 4 fm. In 
this way we can neglect surface effects on this quantity. The error bars indicates the uncertainty 
related to the statistics of the simulations. We remark that in all the calculations concerning 
the present work in the range of density explored no cluster formation has been detected by our 
numerical procedure which was applied at each time step to check for the existence of aggregation 
processes. Only in the lower limit p = 0.7po and for only part of the time, the formation of 
a big cluster with a mass about 98 % the total one and light particles have been observed. 
From the figure it is clearly seen the depletion of the Pi distribution at small distances due to 
Pauli constraining in CoMD. For the other distributions an enhancement is produced due to 
the attractive effect of the two-body interaction. In particular we note that this enhancement 
is more pronounced for P np as compared to Po; this is clearly an effect due to the iso- vectorial 
forces which generate a stronger attraction for neutron-proton couples. Now from eq.(35) we 
can obtain an expression for S v as: 



p = 1 ± koi/(r) 



(33) 




(34) 



P(r) = p(l±k i/(r)) 



(35) 




(36) 
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Figure 1. Probability distribution Pi to find two identical nucleons at a relative distance 
r. In the same figure Pq represents the probability evaluated for neutron and proton couples 
with opposite spin. Finally P np represents the probability obtained for neutron-proton couples. 
(Color) 



1 = 7^1 1 ^° e ^[-z^M^ 

From the above equations we see that even if the overlap integral per nucleon is closely related to 
the density p, it admits correction terms related to the spacial correlation through the integral /. 
In general, for function v(r) localized within a distance of the order of I the value of / decreases 

r 

with the ratio x = ^f: for example, for an exponential profile e~ i the correction terms go to 

-4 

zero like e . In other words, according to what observed in the previous section, when the 
single particles space distribution is enough large an averaging of the spacial correlations effect is 
obtained and the CoMD functional tends to the one represented in eq.(10) which is typical of the 
semi-classical mean-field approximation in transport theory. On the contrary for well localized 
wave-packets with a r of the order of 1-2 fm, I is different from zero and reflects the behavior 
of is(r) at small distances. The above relations can be also generalized for a multi-component 
system characterized by a charge/mass parameter different from zero. 

The determination of the correlation functions v{r) and of the related integral I is a problem 
which can not be solved in a general way. Some special cases are well known. They concern 
the study of the density fluctuations in the hydrodynamical limit valid for large distances 
compared to the mean free-path in non-equilibrium cases (see also |1U| [TT]). However in this 
limit simultaneous spacial correlations are supposed to be zero at small distances. At short 
distances, comparable with the range of the effective interaction, as in our case, approximations 
schemes can be developed actually corresponding to a power decomposition in p as we will 
perform in the next section. 



(37) 



3. The nuclear matter simulation 

As mentioned in the introduction the main goal of the present work is to understand also 
at a quantitative level and for a given simple form of the effective interaction, what are the 
consequences produced by the correlations above discussed by taking as a reference point the 
results obtainable in the framework of the Se-MFA approach. The study will be performed in a 
narrow region of density around the saturation density pq. To this aim in the following we will 
try to find the set of parameters for the effective interaction which reproduce in both the cases 
some of the commonly accepted saturation properties of symmetric nuclear matter. In particular 
we will refer to po=0.165 /to -3 , binding energy E(po) = —16 MeV, compressibility modulus for 
symmetric nuclear mater Knm(po) = 220 MeV. The strength of the iso- vectorial interaction T4 
will be chosen in such way to produce a symmetry energy value €sym (j>o) ~ 28.6 MeV in the 
case of the Se-MFA, i.e. T4 = 32 MeV. From the functional E reported in eqs. (10-12) we can 
find easily the parameter values satisfying the above conditions by solving the system formed 
by the following set of equations for symmetric NM: 

0.165 (38) 
-16MeV (39) 

(40) 
K N m(po) = 220MeV (41) 

In concrete cases due to the finite steps with which we perform the variation on the parameters 
the values of E(po) and Knm(po) are obtained within ±0.5%. Within the above specified 
uncertainty the solution gives the following values for the parameters: To = —263 ±1.3 MeV, 
T 3 = 208 ± 1 MeV and a = 1.25 ± 0.02. 
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3.1. NM calculations and CoMD model 

The evaluation of the total energy per nucleon E c related to the CoMD calculations requires 
the solution of the many-body problem using the equation of motion regulating the wave-packet 
dynamics. Details on the numerical procedure for the CoMD calculations are described in 
ref . [13 E3| - 

At the different densities changing between 0.7-1.2/30 with steps equal to 0.1po> the 
calculations have been performed by enclosing a relatively large number of particles A\ = 1600 
and A2 = 3560 in spherical volumes of radii R = tqA 1 ^ 3 with ro = (j^^) 1 ^ 3 - Particles trying to 
escape from the spheres have been re-scattered inside through an elastic reflection at the surface. 
For symmetrical configuration, starting from the parameter values minimizing the functional E 
in eq. (10-12) we have searched for the stationary minimum energy conditions by applying the 
cooling-warming procedure coupled with the constraint related to the Pauli principle as described 
in Ref.|13j. Calculations have been performed for the two systems having number of particles 
equal to A± and A%. The value of T4 has been fixed to 32 MeV and the calculations have 
been performed for different values of the stiffness parameter 7. From the minimum energy 
configurations we have evaluated the related intensive quantities S v (p), a(p),pA(p) and E^ n . 
Corrections due to the surface effects, which are necessary to estimate the associate bulk values 
have been evaluated using the following relation: 

Qi = Q b + Q s a; 1/3 + oa5Q s a; 2/3 (42) 

Qi indicates the quantity valuated for the system with mass Ai (i=l,2). Qt, and Q s are the 
bulk and surface coefficients. Effects related to the curvature are represented by the last term 




Figure 2. Typical result for the primary quantities S v (p), a(p),pA(p) as a function of reduced 
density. The star symbols indicate the results of the study performed on the lighter system with 
mass Ax, the squares represent the results obtained for heavier system containing A2 particles. 
Finally the circles indicate the obtained corrected values for surface effect according to eq.(42). 
The black solid lines are the final results of a fit procedure with a second order polynomial of 
the density . The red line represents the function p as a function of the reduced density. 



of eq.(42). The coefficient 0.45 has been deduced performing a couple of calculations in boxes 
having the same volume as the considered spheres. 

As an example, for (3 = and 7 = 1 (form factor F2) in Fig. 2 we show as a function of the 
reduced density the values of S v (p), a(p), pa{p) evaluated for the systems 1 (marked points 
with star), for the system 2 (marked points with square) and the bulk estimated values (dot 
points). Corrections of the same order are evaluated also for the kinetic energy contribution. 
In the following for simplicity we will refer to the bulk quantities without using the subscript 
b. After this first step, we fit the evaluated bulk quantities with a polynomial function of the 
reduced density — . The above quantities show in fact deviations from a simple linear behavior 
as a function of p. This can be seen by looking at the red line in the figure {p as a function of 
p/po)- A second order polynomial reproduces very well the behavior in the range of explored 
densities. The results of the fit are shown in the figures with lines. 

The obtained functions are then substituted in E c and therefore the total binding energy can 
be now evaluated with continuity as a function of p and we can search for the parameter values 
solving the CoMD functional E c obtained in Sec. II. 1 and satisfying the conditions expressed 
in eqs. (38-41). We have to note that the numerical solution of this system of coupled equations 
in general can not be obtained with the same precision as the one involving the functional E. 
In most of the cases, depending on the stiffness parameter 7 it is possible to obtain solutions 
reproducing the E(pq) and po values within 10% while a larger spread is obtained for the 
Knm{po)- The chosen solutions will be the one which minimize the total relative difference 
from the reference values. Having found the best solution for the functional E c , in the sense 
above specified, with the new set of parameter values for To, T3 and a we perform another series 
of microscopic NM simulations on the systems of A±,A2 particles. After having included the 
correction for surface effects, we do the polynomial fit. Then using the new calculated quantities 



0,12 



0,7 0.8 0,9 1 1,1 1.2 1.3 

p/Po 

Figure 3. Final values of S v (p) for symmetric NM as a function of reduced density and for 
different form factors Fk as indicated in the legend. The black solid lines are the result of a fit 
procedure with a second order polynomial of the density. The red line represent the density p 

we solve another time the functional E . This iterative procedure is continued until the values 
of the parameters differs, in two subsequent steps, by an amount less than ±5 %. The method 
converges after 2-3 iterations. 

4. Results and discussion 

In this section we illustrate the results obtained from the recursive procedure previously 
described. 

4-1- Results on the primary quantities 

As an example in Fig. 3 we show the final values of S v as a function of the reduced density 
obtained in the case of symmetric NM and for different form factors Fk (T4 = 32 MeV). The 
lines represent the results of the fit with a second order polynomial. The red line represents 
instead the linear relation corresponding to the density p. As we can see S v shows deviations 
from p depending also on the used form factor. 

Under the same conditions in Fig. 4 with black line and points, we show the value of a as a 
function of the reduced density. The red line (and points) represents the values of a obtained in 
the case of T4 = MeV using the form factor F<i- The blue line represents the values obtained 
for T4 = 59 MeV which corresponds in the case of the Se-MFA to a value of e sym (po) equal to 
42 MeV. 

As can be seen the value of a decreases with the density and increases with T4. For T4 = 
finite values of a are essentially due to correlations imposed by the Pauli principle in the system 
of interacting particles. 

4-2. Total binding energy 

In Fig. 5, upper panel, we show with a solid line the total energy per nucleon E as a function of the 
reduced density (eqs.(ll,12)) satisfying the requested conditions on NM saturation properties. 
This result represent our reference point. 
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Figure 4. Final values of the a correlation coefficient for j3 = as a function of reduced 
density. The black points and lines (results of a second order polynomial fit) refers to different 
form factors Fk and T4 = 32 MeV. For the forma factor F% the red and blu points refers to T4 
values equal to and 59 MeV respectively. 

The curves with broken lines represent the results of our NM simulations with CoMD model at 
the first step of the iteration for the different indicated form factors Fk ■ In this case the parameter 
used for the effective interaction are the same like the ones obtained from the minimization of 
the energy functional related to the Se-MFA. As can be seen the curves in this cases are rather 
different. In particular for F\ and F2 the minimum energies are shifted to lower values and 
the compressibility modulus also shows large deviations compared to the chosen reference value. 
F3 shows instead even a maximum at the saturation density. In the lower panel we show the 
corresponding values of Euas (see eq.(32)). Together with the density dependence of the primary 
quantities above described, this term is the main co-responsible of the observed deviations. The 
figure shows that the density behavior of Eu as strongly depends at a quantitative and qualitative 
level on the used form factors Fk- In particular we note an increasing slope from positive to 
negative values with the increasing of the degree of stiffness 7 characterizing the behavior of the 
iso-vectorial forces. We can expect therefore that the necessary correction on the parameters of 
the Iso-scalar effective interaction to reproduce the reference properties of the symmetric NM 
will show a dependence on 7. 

In Fig. 6 we show analogues results after that the self consistent iteration procedure has been 
completed as described in Sec. Ill . In the interval 7 ~ 0.85 — 1.5 the saturation density and 
binding energies are reproduced within some percent. The compressibility is instead obtained 
within about 20%. In the figure we show some example of these solutions with black lines. 
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Figure 5. Upper panel:total energy per nucleon E for symmetric NM as a function of the 
reduced density. The solid line represents the results obtained in the Se-MFA. With discontinue 
lines we plot the results for CoMD calculations corresponding to the first step of the iterative 
procedure (see the text). Different discontinue lines represents results related to different form 
factors -Ffc according to the legend. Bottom panel: For the same parameter values the values of 
Euas (see eq.(32)) are plotted as function of the reduced density. 
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Figure 6. Total Energy per nucleon E obtained through CoMD calculations for symmetric 
NM as a function of the reduced density. Different curves refer to different form factors F^. In 
the legend the value of Knm at Po are also indicated. The red and magenta lines represent an 
example of limiting cases reached when 7 is well outside the range 1.5-0.85 



In general we note a more pronounced asymmetry of the curves around the saturation density 
with a reduced slope of the lower density branch and an increased slope for reduced density 
larger than 1. This behaviors is a direct consequence of the density dependence of the average 
overlap integral (see Fig. 5). For stiffness parameter values out of the indicated interval we 
observe a fast increasing of the compressibility (well beyond 400 MeV) and of the corresponding 
binding energy at the saturation density. In Fig. 6 an example showing this trend is represented 
by colored lines obtained for 7 = 0.75 and 7 = 1.75. All these circumstances shows that, in 
the frame work of the present molecular dynamics approach, the F4 form factor with stiffness 
parameter values external to the range 7 ~ 1.5 — 0.85, due to the correlations discussed, is not 
able to reproduce the chosen commonly accepted saturation properties of the symmetric nuclear 
matter. These results are consistent with recent findings on the study of the 40 ' 48 Ca + 40 > 48 Ca 
systems at 25 MeV/nucleon [22 E2] concerning the balance between the yields of incomplete- 
fusion and multi-breakup processes. Finally in Fig. 7 we show as a function of the 7 parameter 
the values of the To, T3 and a obtained from the iterative procedure described in Sec. III. 1 . 
Apart from the extremal plotted values, corresponding to high values of the compressibility ( 
see Fig.6), the internal ones correspond to saturation density, binding energies and K^m values 
well within 20% of the value obtained in the case of the Se-MFA. From the figure we observe a 
dependence of the parameter values describing the iso-scalar forces on the stiffness parameter 7 
associated to the iso-vectorial interaction. In the internal region of the explored interval, even 
if the dependence is moderate, maximum changes of the order of 16% and 20% are obtained for 
the T3 and a values respectively. However it is remarkable that the average values of the iso- 
scalar interaction parameters show large differences compared to the reference values obtained 
in the case of the Se-MFA (T = -263 MeV,T 3 = 208 MeV and a = 1.25. see Sec.III) which are 
instead independent on 7. 

5. Summary and Concluding Remarks 

Many-body correlations produced in Molecular Dynamics approach based on the CoMD model 
have been discussed and their connection with the used effective interaction have been analyzed 
in the case of symmetric nuclear matter simulations. This study has been performed by 
comparing the results obtained for the total energy functional and the nuclear matter (NM) 
main saturation properties with the ones obtainable in the case of a semiclassical mean field 
approximations (Se-MFA). This comparison has been performed by using the same kind of 
simple effective Skyrme interaction. While we can expect the two approaches to produce large 
differences at low density due to the cluster formation process, the following study shows that 
noticeable differences are obtained also in a narrow range of densities around the saturation 
one where no cluster production has been observed. The effective Skyrme interaction include 
two-body and three-body effective iso-scalar interactions plus a two-body iso-vectorial one with 
form factors commonly used also in Se-MFA. In CoMD model calculations, around the saturation 
density, the effects related to the spacial correlations generated by the usage of the localized wave- 
packets and the ones associated to the multi-particles correlations produced through the Pauli 
principle constraint can be well described through a second order polynomial decomposition of 
the total energy as a function of the density. The obtained results show that, contrary to the case 
of the Se-MFA, the discussed correlations produce an interdependence between the parameters 
describing the iso-scalar forces and the ones related to the iso-vectorial interaction. The usage of 
an iterative procedure tuned to obtain in both the cases (Se-MFA and in CoMD approaches) very 
similar saturation density, total energy and compressibility values for symmetric NM allowed 
to extract the "good" set of parameters used for CoMD calculations. The obtained values 
differ under many aspects from the ones obtained from the Se-MFA: in particular the density 
dependence of the used form factors describing the iso-vectorial forces has to change now in a 
more restricted range of stiffness values. Moreover the values of the coefficients describing the 




Figure 7. The values of the parameters To T3 and a (panel a,b,c respectively) obtained through 
the iterative procedure applied to CoMD calculations (see the Sec. III.l)are shown as a function 
of the stiffness parameter 7. The values of 7 equal to 1 and 1.5 are associated to the functional 
form F2 and F\ respectively (see Sec. II). The others ones are instead associated to F4. The 
bar errors represent the global uncertainty on the parameter values related to the numerical 
procedures 



iso-scalar interactions are rather different in the two cases. Work is in progress to extend these 
studies to asymmetric NM. Finally we conclude by observing that even if from a numerical point 
of view the obtained results are strictly valid for the CoMD model, the performed study shows 
that the observed differences in the parameter values describing the chosen effective interaction 
can have a more wide meaning. They are indeed strictly linked to some general properties of 
the semiclassical wave packets dynamics. 
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